{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Levelset Segmentation <a href=\"https://mybinder.org/v2/gh/InsightSoftwareConsortium/SimpleITK-Notebooks/master?filepath=Python%2F31_Levelset_Segmentation.ipynb\"><img style=\"float: right;\" src=\"https://mybinder.org/badge_logo.svg\"></a>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import SimpleITK as sitk\n",
    "from myshow import myshow, myshow3d\n",
    "\n",
    "# Download data to work on\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_T1 = sitk.ReadImage(fdata(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd\"))\n",
    "img_T2 = sitk.ReadImage(fdata(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd\"))\n",
    "img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)\n",
    "img_T1_255 = sitk.Cast(sitk.RescaleIntensity(img_T1), sitk.sitkUInt8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "sitk.Show(img_T1, title=\"T1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = (106, 116, 141)\n",
    "pt = img_T1.TransformIndexToPhysicalPoint(idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seg = sitk.Image(img_T1.GetSize(), sitk.sitkUInt8)\n",
    "seg.CopyInformation(img_T1)\n",
    "seg[idx] = 1\n",
    "seg = sitk.BinaryDilate(seg, [3] * 3)\n",
    "myshow3d(\n",
    "    sitk.LabelOverlay(img_T1_255, seg),\n",
    "    zslices=range(idx[2] - 3, idx[2] + 4, 3),\n",
    "    dpi=30,\n",
    "    title=\"Initial Seed\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = sitk.LabelStatisticsImageFilter()\n",
    "stats.Execute(img_T1, seg)\n",
    "print(stats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factor = 1.5\n",
    "lower_threshold = stats.GetMean(1) - factor * stats.GetSigma(1)\n",
    "upper_threshold = stats.GetMean(1) + factor * stats.GetSigma(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_ls = sitk.SignedMaurerDistanceMap(seg, insideIsPositive=True, useImageSpacing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsFilter = sitk.ThresholdSegmentationLevelSetImageFilter()\n",
    "lsFilter.SetLowerThreshold(lower_threshold)\n",
    "lsFilter.SetUpperThreshold(upper_threshold)\n",
    "lsFilter.SetMaximumRMSError(0.02)\n",
    "lsFilter.SetNumberOfIterations(100)\n",
    "lsFilter.SetCurvatureScaling(1)\n",
    "lsFilter.SetPropagationScaling(1)\n",
    "lsFilter.ReverseExpansionDirectionOn()\n",
    "ls = lsFilter.Execute(init_ls, sitk.Cast(img_T1, sitk.sitkFloat32))\n",
    "print(lsFilter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zslice_offset = 4\n",
    "t = \"LevelSet after \" + str(lsFilter.GetNumberOfIterations()) + \" iterations\"\n",
    "myshow3d(\n",
    "    sitk.LabelOverlay(img_T1_255, ls > 0),\n",
    "    zslices=range(idx[2] - zslice_offset, idx[2] + zslice_offset + 1, zslice_offset),\n",
    "    dpi=20,\n",
    "    title=t,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsFilter.SetNumberOfIterations(25)\n",
    "img_T1f = sitk.Cast(img_T1, sitk.sitkFloat32)\n",
    "ls = init_ls\n",
    "niter = 0\n",
    "for i in range(0, 10):\n",
    "    ls = lsFilter.Execute(ls, img_T1f)\n",
    "    niter += lsFilter.GetNumberOfIterations()\n",
    "    t = (\n",
    "        \"LevelSet after \"\n",
    "        + str(niter)\n",
    "        + \" iterations and RMS \"\n",
    "        + str(lsFilter.GetRMSChange())\n",
    "    )\n",
    "    fig = myshow3d(\n",
    "        sitk.LabelOverlay(img_T1_255, ls > 0),\n",
    "        zslices=range(\n",
    "            idx[2] - zslice_offset, idx[2] + zslice_offset + 1, zslice_offset\n",
    "        ),\n",
    "        dpi=20,\n",
    "        title=t,\n",
    "    )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
